import numpy as np
from fractions import Fraction
import math
import matplotlib.pyplot as plt


def gauss_elimination(matrix, vector):
    if len(matrix) != len(vector):
        print('ERROR: the size of matrix and vector does not match!')
        return
    M = np.column_stack([matrix, vector])
    n = len(matrix)
    for i in range(n-1):
        for j in range(i+1, n):
            factor = Fraction(M[j, i], M[i, i])
            M[j, i:] -= factor * M[i, i:]
    x = np.zeros(n, Fraction)
    for i in range(n-1, -1, -1):
        x[i] = (M[i, n] - np.dot(M[i, i:n], x[i:n])) / M[i, i]
    return x


def hilbert_matrix(n):
    H = np.zeros((n, n), Fraction)
    for i in range(n):
        for j in range(n):
            H[i][j] = Fraction(1, i+j+1)
    return H


def get_r_parted(n):
    c1 = np.zeros(n, Fraction)
    c2 = np.zeros(n, Fraction)
    for i in range(n):
        c1[i] = pow(-1, i)*Fraction(1, 1)
        c2[i] = Fraction(0, 1)
        for j in range(i):
            c2[i] += pow(-1, i+j+1)*Fraction(1, j+1)
    return c1, c2


def get_beta_gama(n):
    Hilbert_n = hilbert_matrix(n)
    c1, c2 = get_r_parted(n)
    beta = gauss_elimination(Hilbert_n, c1)
    gama = gauss_elimination(Hilbert_n, c2)
    print('n = ', n, ': ')
    print('cond(Hn) = ', np.linalg.cond(Hilbert_n.astype(float)))
    print('beta:', end=' ')
    for i in range(n):
        print(beta[i].numerator, '/', beta[i].denominator, end='    ')
    print()
    print('gama:', end=' ')
    for i in range(n):
        print(gama[i].numerator, '/', gama[i].denominator, end='    ')
    print()
    return beta, gama


def output(matrix):
    s = [[str(e) for e in row] for row in matrix]
    lens = [max(map(len, col)) for col in zip(*s)]
    fmt = "\t".join("{{:{}}}".format(x) for x in lens)
    table = [fmt.format(*row) for row in s]
    print("\n".join(table))


ln2 = 0.69315
Tik_para = 1e-10

alpha_mp = np.zeros((6, 7))
alpha_fp = np.zeros((6, 7))
alpha_tk = np.zeros((6, 7))
error_mp = np.zeros((6, 7))
error_fp = np.zeros((6, 7))
error_tk = np.zeros((6, 7))

for n in range(1, 7):
    H_n_ra = hilbert_matrix(n+1)
    H_n_fl = H_n_ra.astype(float)
    H_n_tk = Tik_para*np.eye(n+1) + np.dot(np.transpose(H_n_fl), H_n_fl)

    r_n_part1, r_n_part2 = get_r_parted(n+1)
    r_n = r_n_part1*ln2 + r_n_part2
    r_n_fl = r_n.astype(float)
    r_n_tk = np.dot(np.transpose(H_n_fl), r_n_fl)

    beta, gama = get_beta_gama(n+1)
    print()

    alpha_mp[n-1, 0:n+1] = math.log(2)*beta + gama
    alpha_fp[n-1, 0:n+1] = np.linalg.solve(H_n_fl, r_n_fl)
    alpha_tk[n-1, 0:n+1] = np.linalg.solve(H_n_tk, r_n_tk)

    for i in range(n+1):
        error_mp[n-1, i] = abs((pow(-1, i)-alpha_mp[n-1, i])/pow(-1, i))
        error_fp[n-1, i] = abs((pow(-1, i)-alpha_fp[n-1, i])/pow(-1, i))
        error_tk[n-1, i] = abs((pow(-1, i)-alpha_tk[n-1, i])/pow(-1, i))

print('(3) formula computation by approximating ln2 to machine precision: ')
output(alpha_mp)
print()
print('(4) matrix calculation by floating-point arithmetic with ln 2 ≈ 0.69315: ')
output(alpha_fp)
print()
print('(5) matrix calculation by Tikhonov regularization with α = 1e−10: ')
output(alpha_tk)
print()

x_axis = [1, 2, 3, 4, 5, 6]
for i in range(3):
    plt.xlabel('polynomial order')
    plt.ylabel('error with Taylor expansion coefficients')
    plt.plot(x_axis, error_mp[:, i], color='b')
    plt.plot(x_axis, error_fp[:, i], color='g')
    plt.plot(x_axis, error_tk[:, i], color='y')
    plt.legend(('formula computation',
                'matrix calculation (floating-point arithmetic)',
                'matrix calculation (Tikhonov regularization)'), loc='upper left')
    plt.savefig('./pic/convergence_of_a'+str(i))
    plt.show()
